Read in the packages. The working directory is wherever the R Notebook is located.
Read in the spreadsheet and take a look at the data.
###read in spreadsheet
loc <- read_xlsx("Original-spreadsheets/all species New_6-14-19.xlsx") %>%
janitor::clean_names() %>%
mutate(reproductive_mode = as.factor(reproductive_mode))
#get the number of individuals, and the sexuality counts per species
count_repro_mode <- loc %>%
group_by(genus, species, reproductive_mode) %>%
count() %>%
mutate(genus_species = str_c(genus, species, sep = "_"),
genus_species = str_replace_all(genus_species, " ", "_"),
genus_species = str_replace_all(genus_species, "\\.", "")) %>%
ungroup() %>%
mutate(genus_species = fct_reorder(genus_species, n, sum)) %>%
ggplot(aes(x = genus_species, y = n, fill = reproductive_mode)) +
geom_col() +
coord_flip() +
theme_minimal()
count_repro_mode##Map Plot a leaflet map of the localities. The leaflet map is interactive. You can click on the localities and a flag with some metadata will pop up!
#make locality shape file and assign WGS coord system
coord_points <- st_as_sf(loc, coords = c("longitude", "latitude"),
crs = 4326, agr = "constant")
#use sourced plot_locs_leaflet script to plot localities
all_plot <- plot_locs_leaflet(loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(all_plot, url = paste0(getwd(), "/plots/repro_mode_plots/all_species_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/all_species_map.pdf"))##PCA-Genera {.tabset}
###Climate Data Obtain the bioclim layers for analysis. I’m using all 19 for this preliminary exploration. I plotted the first bioclim just to make sure nothing seems wonky. I’m using CHELSA data downloaded from their website. Since the files are huge, I only unzip them one at a time, crop them, and write them to GeoTiff files that I can then load in as a rasterstack.
##get chelsa data
#chelsa_folder <- "/Users/connorfrench/Dropbox/Old_Mac/climate-data/chelsa_30s_bio"
#zip_files <- list.files(chelsa_folder, full.names = TRUE)
#using the Unarchiver commandline tools for Mac to unzip the 7zip chelsa layers. Regular unzip() does not work with 7z zipped files
#for (file in zip_files) {
#set temp directory
# tempd <- tempdir()
# system(paste("unar", file, "-o", tempd))
# r <- raster(list.files(tempd, pattern = "*.tif", full.names = TRUE)) %>%
# crop(extent(166, 179, -48, -34))
# writeRaster(r, filename = paste0("~/Desktop/", list.files(tempd, pattern = "*.tif")), format = "GTiff")
# unlink(tempd, recursive = TRUE)
#}
clim_files <- "/Users/connorfrench/Dropbox/Old_Mac/climate-data/chelsa_30sec_NewZealand/chelsa_bioclims_NZ.tif"
w <- stack(clim_files)###PCA by locality This is a PCA of the climate data extracted for each locality, rather than a PCA of the total climate space.
Run the pca and check out variable loadings and proportion of variance explained by components.
#extract data from worldclim for each locality. Making this into a data frame with columns labeled so the row labeling lines up after I remove the NAs.
#extract data from worldclim for each locality.
coords <- data.frame(latitude = loc$longitude, longitude = loc$latitude)
loc.clim <- dplyr::bind_cols(loc, raster::extract(w, coords, method = "simple", df = TRUE)) %>%
drop_na(chelsa_bioclims_NZ.1) %>%
dplyr::select(-ID)
#make a matrix of only bioclim values
clim.mat <- loc.clim[,grep("bio", names(loc.clim))] %>% as.matrix()
#run pca on climate variables
clim.pca <- prcomp(clim.mat, scale = TRUE)
summary.pca <- summary(clim.pca) #check out the components
#plot tables
summary.pcaImportance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
Standard deviation 2.9511 2.4947 1.36934 1.02403 0.80149 0.54359 0.32182 0.24972 0.11677 0.10148
Proportion of Variance 0.4584 0.3276 0.09869 0.05519 0.03381 0.01555 0.00545 0.00328 0.00072 0.00054
Cumulative Proportion 0.4584 0.7859 0.88462 0.93981 0.97362 0.98917 0.99462 0.99790 0.99862 0.99916
PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.08340 0.05671 0.04746 0.03910 0.03147 0.02392 0.01408 0.01153 0.00911
Proportion of Variance 0.00037 0.00017 0.00012 0.00008 0.00005 0.00003 0.00001 0.00001 0.00000
Cumulative Proportion 0.99953 0.99970 0.99982 0.99990 0.99995 0.99998 0.99999 1.00000 1.00000
| PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | -0.301 | 0.148 | 0.174 |
| chelsa_bioclims_NZ.2 | -0.285 | 0.143 | 0.276 |
| chelsa_bioclims_NZ.3 | -0.310 | 0.142 | 0.097 |
| chelsa_bioclims_NZ.4 | -0.150 | -0.358 | 0.010 |
| chelsa_bioclims_NZ.5 | -0.189 | -0.327 | 0.024 |
| chelsa_bioclims_NZ.6 | -0.135 | -0.362 | -0.003 |
| chelsa_bioclims_NZ.7 | -0.169 | 0.136 | 0.089 |
| chelsa_bioclims_NZ.8 | -0.188 | -0.328 | 0.031 |
| chelsa_bioclims_NZ.9 | -0.130 | -0.365 | -0.006 |
| chelsa_bioclims_NZ.10 | -0.119 | -0.361 | -0.036 |
| chelsa_bioclims_NZ.11 | -0.214 | -0.280 | 0.094 |
| chelsa_bioclims_NZ.12 | 0.277 | -0.106 | 0.360 |
| chelsa_bioclims_NZ.13 | 0.261 | -0.105 | 0.306 |
| chelsa_bioclims_NZ.14 | 0.248 | -0.083 | 0.383 |
| chelsa_bioclims_NZ.15 | -0.219 | 0.123 | 0.492 |
| chelsa_bioclims_NZ.16 | -0.314 | 0.143 | 0.009 |
| chelsa_bioclims_NZ.17 | 0.272 | -0.106 | 0.380 |
| chelsa_bioclims_NZ.18 | -0.167 | 0.128 | -0.084 |
| chelsa_bioclims_NZ.19 | -0.247 | 0.100 | 0.317 |
Two plots: One plot of the PCA colored according to genus, with convex hulls surrounding the genera. It looks like this reflects a latitudinal gradient in temperature! You can interact with the PCA plot by clicking on points to view associated metadata. You can isolate the genus you want to view by double clicking the genus in the legend! You can also remove a genus by clicking on it once. There’s some other functionality you can explore in the toolbar at the top of the plot. The second plot is a PCA colored according to reproductive mode. It looks like asexual populations occupy slightly larger niche space, but both reproductive modes have a similar niche center.
#add pca results to loc.clim data frame
loc.clim <- data.frame(loc.clim, clim.pca$x)
#use sourced plot_clim_pca function to plot the pca results. args are the data set with species names and PC axis values and the pca summary
all_pca <- plot_clim_pca(loc.clim, summary.pca, factor = "genus")Ignoring unknown aesthetics: text
#use sourced plot_clim_pca function to plot the pca results. args are the data set with species names and PC axis values and the pca summary
repro_pca <- plot_clim_pca(loc.clim, summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#save the plot colored by genus
#htmlwidgets::saveWidget(all_pca, paste0(getwd(), "/plots/repro_mode_plots/all_species_pca_genus.html"), selfcontained = TRUE)
#save the plot colored by reproductive mode
#htmlwidgets::saveWidget(repro_pca, paste0(getwd(), "/plots/repro_mode_plots/all_species_pca_repro.html"), selfcontained = TRUE)##PCA-Species {.tabset} These are PCAs of environmental space for species within genera. Each climate PCA is of localities for a single genus, colored by species. I’m doing this even for genera with one species, so it’s easy to see if certain localities group together.
###Acanthoxyla
#source function to conduct a PCA on individual species
summary.list.acan <- species_pca_fun(loc.clim, "acanthoxyla")
#plot
acan_plot <- plot_clim_pca(summary.list.acan$loc.clim, summary.list.acan$summary.pca, "reproductive_mode")Ignoring unknown aesthetics: text
#save pca plot
#htmlwidgets::saveWidget(acan_plot, paste0(getwd(), "/plots/repro_mode_plots/acanthoxyla_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
acan_loc <- loc %>%
filter(genus == "acanthoxyla")
#use sourced plot_locs_leaflet script to plot localities
acan_map <- plot_locs_leaflet(acan_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(acan_map, url = paste0(getwd(), "/plots/repro_mode_plots/acan_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/acan_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
Standard deviation 2.8059 2.5135 1.7771 0.88976 0.66876 0.44564 0.36423 0.22139 0.1064 0.09506
Proportion of Variance 0.4144 0.3325 0.1662 0.04167 0.02354 0.01045 0.00698 0.00258 0.0006 0.00048
Cumulative Proportion 0.4144 0.7469 0.9131 0.95475 0.97829 0.98874 0.99573 0.99831 0.9989 0.99938
PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.06761 0.05169 0.03848 0.03336 0.02780 0.02371 0.01891 0.01400 0.01128
Proportion of Variance 0.00024 0.00014 0.00008 0.00006 0.00004 0.00003 0.00002 0.00001 0.00001
Cumulative Proportion 0.99962 0.99976 0.99984 0.99989 0.99993 0.99996 0.99998 0.99999 1.00000
loadings.acan <- summary.list.acan$summary.pca$rotation
knitr::kable(round(loadings.acan[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | -0.266 | -0.230 | 0.170 |
| chelsa_bioclims_NZ.2 | -0.238 | -0.200 | 0.297 |
| chelsa_bioclims_NZ.3 | -0.277 | -0.240 | 0.056 |
| chelsa_bioclims_NZ.4 | 0.262 | -0.267 | 0.039 |
| chelsa_bioclims_NZ.5 | 0.236 | -0.288 | 0.068 |
| chelsa_bioclims_NZ.6 | 0.268 | -0.258 | 0.015 |
| chelsa_bioclims_NZ.7 | -0.192 | -0.051 | 0.217 |
| chelsa_bioclims_NZ.8 | 0.239 | -0.286 | 0.073 |
| chelsa_bioclims_NZ.9 | 0.269 | -0.256 | 0.012 |
| chelsa_bioclims_NZ.10 | 0.272 | -0.249 | -0.011 |
| chelsa_bioclims_NZ.11 | 0.191 | -0.283 | 0.148 |
| chelsa_bioclims_NZ.12 | 0.219 | 0.232 | 0.285 |
| chelsa_bioclims_NZ.13 | 0.232 | 0.246 | 0.151 |
| chelsa_bioclims_NZ.14 | 0.128 | 0.138 | 0.463 |
| chelsa_bioclims_NZ.15 | -0.150 | -0.097 | 0.477 |
| chelsa_bioclims_NZ.16 | -0.278 | -0.245 | -0.034 |
| chelsa_bioclims_NZ.17 | 0.197 | 0.197 | 0.370 |
| chelsa_bioclims_NZ.18 | -0.074 | -0.304 | -0.109 |
| chelsa_bioclims_NZ.19 | -0.246 | -0.049 | 0.321 |
###Argosarchus
#conduct pca
summary.list.argo <- species_pca_fun(loc.clim, "argosarchus")
#plot
argo_plot <- plot_clim_pca(summary.list.argo$loc.clim, summary.list.argo$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(argo_plot, paste0(getwd(), "/plots/repro_mode_plots/argosarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
argo_loc <- loc %>%
filter(genus == "argosarchus")
#use sourced plot_locs_leaflet script to plot localities
argo_map <- plot_locs_leaflet(argo_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(argo_map, url = paste0(getwd(), "/plots/repro_mode_plots/argo_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/argo_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
Standard deviation 2.8508 2.4341 1.6389 1.09192 0.76813 0.51507 0.35817 0.2222 0.10883 0.10035
Proportion of Variance 0.4277 0.3118 0.1414 0.06275 0.03105 0.01396 0.00675 0.0026 0.00062 0.00053
Cumulative Proportion 0.4277 0.7396 0.8809 0.94367 0.97473 0.98869 0.99544 0.9980 0.99866 0.99919
PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.08126 0.05441 0.04318 0.03895 0.03285 0.02417 0.01722 0.01643 0.01134
Proportion of Variance 0.00035 0.00016 0.00010 0.00008 0.00006 0.00003 0.00002 0.00001 0.00001
Cumulative Proportion 0.99954 0.99970 0.99980 0.99988 0.99993 0.99996 0.99998 0.99999 1.00000
loadings.argo <- summary.list.argo$summary.pca$rotation
knitr::kable(round(loadings.argo[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | -0.266 | -0.227 | -0.182 |
| chelsa_bioclims_NZ.2 | -0.241 | -0.208 | -0.293 |
| chelsa_bioclims_NZ.3 | -0.279 | -0.228 | -0.110 |
| chelsa_bioclims_NZ.4 | -0.225 | 0.312 | -0.037 |
| chelsa_bioclims_NZ.5 | -0.257 | 0.272 | -0.023 |
| chelsa_bioclims_NZ.6 | -0.202 | 0.330 | -0.026 |
| chelsa_bioclims_NZ.7 | -0.141 | -0.175 | 0.021 |
| chelsa_bioclims_NZ.8 | -0.255 | 0.275 | -0.033 |
| chelsa_bioclims_NZ.9 | -0.197 | 0.334 | -0.030 |
| chelsa_bioclims_NZ.10 | -0.176 | 0.347 | 0.013 |
| chelsa_bioclims_NZ.11 | -0.277 | 0.219 | -0.090 |
| chelsa_bioclims_NZ.12 | 0.266 | 0.149 | -0.318 |
| chelsa_bioclims_NZ.13 | 0.267 | 0.108 | -0.280 |
| chelsa_bioclims_NZ.14 | 0.227 | 0.162 | -0.310 |
| chelsa_bioclims_NZ.15 | -0.113 | -0.143 | -0.516 |
| chelsa_bioclims_NZ.16 | -0.291 | -0.224 | -0.009 |
| chelsa_bioclims_NZ.17 | 0.258 | 0.163 | -0.324 |
| chelsa_bioclims_NZ.18 | -0.098 | -0.073 | 0.247 |
| chelsa_bioclims_NZ.19 | -0.191 | -0.155 | -0.384 |
Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
###Asteliaphasma
#pca
summary.list.aste <- species_pca_fun(loc.clim, "asteliaphasma")
#plot
aste_plot <- plot_clim_pca(summary.list.aste$loc.clim, summary.list.aste$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(aste_plot, paste0(getwd(), "/plots/repro_mode_plots/asteliaphasma_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
aste_loc <- loc %>%
filter(genus == "asteliaphasma")
#use sourced plot_locs_leaflet script to plot localities
aste_map <- plot_locs_leaflet(aste_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(aste_map, url = paste0(getwd(), "/plots/repro_mode_plots/aste_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/aste_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
Standard deviation 3.104 2.4714 1.27720 0.87566 0.76781 0.38573 0.19883 0.16377 0.1382 0.10580
Proportion of Variance 0.507 0.3215 0.08586 0.04036 0.03103 0.00783 0.00208 0.00141 0.0010 0.00059
Cumulative Proportion 0.507 0.8284 0.91427 0.95463 0.98566 0.99349 0.99557 0.99698 0.9980 0.99858
PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.09487 0.09187 0.06475 0.04986 0.03308 0.03192 0.02037 0.01726 0.01092
Proportion of Variance 0.00047 0.00044 0.00022 0.00013 0.00006 0.00005 0.00002 0.00002 0.00001
Cumulative Proportion 0.99905 0.99949 0.99971 0.99984 0.99990 0.99996 0.99998 0.99999 1.00000
loadings.aste <- summary.list.aste$summary.pca$rotation
knitr::kable(round(loadings.aste[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | -0.298 | -0.090 | -0.229 |
| chelsa_bioclims_NZ.2 | -0.288 | -0.043 | -0.320 |
| chelsa_bioclims_NZ.3 | -0.298 | -0.118 | -0.173 |
| chelsa_bioclims_NZ.4 | 0.198 | -0.312 | -0.110 |
| chelsa_bioclims_NZ.5 | 0.134 | -0.357 | -0.040 |
| chelsa_bioclims_NZ.6 | 0.209 | -0.281 | -0.137 |
| chelsa_bioclims_NZ.7 | -0.167 | -0.205 | 0.152 |
| chelsa_bioclims_NZ.8 | 0.143 | -0.352 | -0.071 |
| chelsa_bioclims_NZ.9 | 0.214 | -0.276 | -0.148 |
| chelsa_bioclims_NZ.10 | 0.187 | -0.284 | -0.161 |
| chelsa_bioclims_NZ.11 | 0.164 | -0.334 | -0.085 |
| chelsa_bioclims_NZ.12 | 0.250 | 0.211 | -0.238 |
| chelsa_bioclims_NZ.13 | 0.230 | 0.192 | -0.354 |
| chelsa_bioclims_NZ.14 | 0.266 | 0.202 | -0.053 |
| chelsa_bioclims_NZ.15 | -0.103 | 0.170 | -0.654 |
| chelsa_bioclims_NZ.16 | -0.295 | -0.154 | -0.044 |
| chelsa_bioclims_NZ.17 | 0.257 | 0.216 | -0.192 |
| chelsa_bioclims_NZ.18 | -0.297 | -0.125 | -0.157 |
| chelsa_bioclims_NZ.19 | -0.211 | -0.027 | -0.154 |
Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
###Clitarchus
summary.list.clita <- species_pca_fun(loc.clim, "clitarchus")
clita_plot <- plot_clim_pca(summary.list.clita$loc.clim, summary.list.clita$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(clita_plot, paste0(getwd(), "/plots/repro_mode_plots/clitarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
clita_loc <- loc %>%
filter(genus == "clitarchus")
#use sourced plot_locs_leaflet script to plot localities
clita_map <- plot_locs_leaflet(clita_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(clita_map, url = paste0(getwd(), "/plots/repro_mode_plots/clita_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/clita_map.pdf"))```r
summary.list.clita$summary.pca
loadings.clita <- summary.list.clita$summary.pca$rotation
knitr::kable(round(loadings.clita[,1:3],3)) #Table of loading scores for the first 3 PCs.
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
Now I'm going to to environmental niche factor analysis between sexual and asexual populations within the species.
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuI2dldCBiYWNrZ3JvdW5kIGVudid0IGZvciB0aGUgc3BlY2llc1xuY2hvb19iZ19lbnYgPC0gYmdfZW52X2Nyb3AoY2xpdGFfbG9jLCBcbiAgICAgICAgICAgICAgICAgICAgICAgICAgIHNwZWNpZXMgPSBcImhvb2tlcmlcIixcbiAgICAgICAgICAgICAgICAgICAgICAgICAgIGVudmlyb25tZW50ID0gdywgXG4gICAgICAgICAgICAgICAgICAgICAgICAgICBidWZmZXIgPSAwLjUpXG5cbiNlbmZhIGZvciB0aGUgc2V4dWFsIHNwZWNpZXNcbmNob29fc2V4dWFsX2VuZmEgPC0gZW5mYV9jYWxjX2Z1bihsb2NzID0gY2xpdGFfbG9jLCBcbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzcGVjaWVzID0gXCJob29rZXJpXCIsIFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHJlcHJvZHVjdGl2ZV9tb2RlID0gXCJzZXh1YWxcIiwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWFza19yYXN0ZXIgPSBjaG9vX2JnX2VudilcblxuI2VuZmEgZm9yIHRoZSBhc2V4dWFsIHNwZWNpZXNcbmNob29fYXNleHVhbF9lbmZhIDwtIGVuZmFfY2FsY19mdW4obG9jcyA9IGNsaXRhX2xvYywgXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHNwZWNpZXMgPSBcImhvb2tlcmlcIiwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHJlcHJvZHVjdGl2ZV9tb2RlID0gXCJhc2V4dWFsXCIsIFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtYXNrX3Jhc3RlciA9IGNob29fYmdfZW52KVxuXG5cbiNwbG90IHRoZSBtYXJnaW5hbGl0eSBzY29yZXNcbm1hcmdpbmFsaXR5X2xvbGxpcG9wKHNleF9tYXJnID0gY2hvb19zZXh1YWxfZW5mYSRtLCBcbiAgICAgICAgICAgICAgICAgICAgYXNleF9tYXJnID0gY2hvb19hc2V4dWFsX2VuZmEkbSxcbiAgICAgICAgICAgICAgICAgICAgZnVsbF9zcGVjaWVzX25hbWUgPSBcIkNsaXRhcmNodXMgaG9va2VyaVwiKVxuXG5gYGAifQ== -->
```r
#get background env't for the species
choo_bg_env <- bg_env_crop(clita_loc,
species = "hookeri",
environment = w,
buffer = 0.5)
#enfa for the sexual species
choo_sexual_enfa <- enfa_calc_fun(locs = clita_loc,
species = "hookeri",
reproductive_mode = "sexual",
mask_raster = choo_bg_env)
#enfa for the asexual species
choo_asexual_enfa <- enfa_calc_fun(locs = clita_loc,
species = "hookeri",
reproductive_mode = "asexual",
mask_raster = choo_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = choo_sexual_enfa$m,
asex_marg = choo_asexual_enfa$m,
full_species_name = "Clitarchus hookeri")
###Micrarchus
summary.list.micra <- species_pca_fun(loc.clim, "micrarchus")
micra_plot <- plot_clim_pca(summary.list.micra$loc.clim, summary.list.micra$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(micra_plot, paste0(getwd(), "/plots/repro_mode_plots/micrarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
micra_loc <- loc %>%
filter(genus == "micrarchus")
#use sourced plot_locs_leaflet script to plot localities
micra_map <- plot_locs_leaflet(micra_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
mapview::mapshot(micra_map, url = paste0(getwd(), "/plots/repro_mode_plots/micra_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/micra_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
Standard deviation 3.4431 2.0743 1.31392 0.94676 0.29817 0.2263 0.20567 0.13292 0.09179 0.08909
Proportion of Variance 0.6239 0.2265 0.09086 0.04718 0.00468 0.0027 0.00223 0.00093 0.00044 0.00042
Cumulative Proportion 0.6239 0.8504 0.94126 0.98843 0.99311 0.9958 0.99803 0.99896 0.99941 0.99983
PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.03357 0.02983 0.02512 0.01449 0.01179 0.01068 0.009585 0.008206 0.007186
Proportion of Variance 0.00006 0.00005 0.00003 0.00001 0.00001 0.00001 0.000000 0.000000 0.000000
Cumulative Proportion 0.99988 0.99993 0.99996 0.99998 0.99998 0.99999 0.999990 1.000000 1.000000
loadings.micra <- summary.list.micra$summary.pca$rotation
knitr::kable(round(loadings.micra[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | 0.268 | 0.116 | -0.222 |
| chelsa_bioclims_NZ.2 | 0.266 | 0.086 | -0.270 |
| chelsa_bioclims_NZ.3 | 0.268 | 0.143 | -0.178 |
| chelsa_bioclims_NZ.4 | -0.242 | 0.257 | -0.101 |
| chelsa_bioclims_NZ.5 | -0.237 | 0.254 | -0.107 |
| chelsa_bioclims_NZ.6 | -0.239 | 0.251 | -0.153 |
| chelsa_bioclims_NZ.7 | 0.145 | -0.116 | -0.015 |
| chelsa_bioclims_NZ.8 | -0.237 | 0.254 | -0.118 |
| chelsa_bioclims_NZ.9 | -0.239 | 0.253 | -0.144 |
| chelsa_bioclims_NZ.10 | -0.240 | 0.250 | -0.147 |
| chelsa_bioclims_NZ.11 | -0.229 | 0.252 | -0.102 |
| chelsa_bioclims_NZ.12 | -0.167 | -0.323 | -0.354 |
| chelsa_bioclims_NZ.13 | -0.140 | -0.310 | -0.429 |
| chelsa_bioclims_NZ.14 | -0.195 | -0.325 | -0.202 |
| chelsa_bioclims_NZ.15 | 0.247 | 0.004 | -0.397 |
| chelsa_bioclims_NZ.16 | 0.268 | 0.170 | -0.110 |
| chelsa_bioclims_NZ.17 | -0.190 | -0.311 | -0.294 |
| chelsa_bioclims_NZ.18 | 0.210 | 0.226 | -0.306 |
| chelsa_bioclims_NZ.19 | 0.267 | 0.119 | -0.185 |
###Niveaphasma
summary.list.nive <- species_pca_fun(loc.clim, "niveaphasma")
nive_plot <- plot_clim_pca(summary.list.nive$loc.clim, summary.list.nive$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(nive_plot, paste0(getwd(), "/plots/repro_mode_plots/niveaphasma_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
nive_loc <- loc %>%
filter(genus == "niveaphasma")
#use sourced plot_locs_leaflet script to plot localities
nive_map <- plot_locs_leaflet(nive_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(nive_map, url = paste0(getwd(), "/plots/repro_mode_plots/nive_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/nive_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
Standard deviation 2.9922 2.4493 1.5457 0.94013 0.73205 0.39659 0.24024 0.11345 0.06508 0.04982
Proportion of Variance 0.4712 0.3157 0.1258 0.04652 0.02821 0.00828 0.00304 0.00068 0.00022 0.00013
Cumulative Proportion 0.4712 0.7870 0.9127 0.95923 0.98743 0.99571 0.99875 0.99943 0.99965 0.99978
PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.04382 0.02749 0.02222 0.01856 0.01631 0.01216 0.01109 0.008331 0.006935
Proportion of Variance 0.00010 0.00004 0.00003 0.00002 0.00001 0.00001 0.00001 0.000000 0.000000
Cumulative Proportion 0.99988 0.99992 0.99995 0.99997 0.99998 0.99999 0.99999 1.000000 1.000000
loadings.nive <- summary.list.nive$summary.pca$rotation
knitr::kable(round(loadings.nive[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | 0.091 | 0.374 | 0.184 |
| chelsa_bioclims_NZ.2 | -0.002 | 0.341 | 0.347 |
| chelsa_bioclims_NZ.3 | 0.156 | 0.358 | 0.053 |
| chelsa_bioclims_NZ.4 | 0.297 | -0.154 | 0.159 |
| chelsa_bioclims_NZ.5 | 0.294 | -0.151 | 0.172 |
| chelsa_bioclims_NZ.6 | 0.297 | -0.158 | 0.135 |
| chelsa_bioclims_NZ.7 | -0.196 | 0.021 | 0.110 |
| chelsa_bioclims_NZ.8 | 0.296 | -0.143 | 0.177 |
| chelsa_bioclims_NZ.9 | 0.298 | -0.158 | 0.146 |
| chelsa_bioclims_NZ.10 | 0.300 | -0.143 | 0.132 |
| chelsa_bioclims_NZ.11 | 0.313 | -0.119 | 0.116 |
| chelsa_bioclims_NZ.12 | -0.234 | -0.233 | 0.274 |
| chelsa_bioclims_NZ.13 | -0.216 | -0.203 | 0.338 |
| chelsa_bioclims_NZ.14 | -0.241 | -0.238 | 0.233 |
| chelsa_bioclims_NZ.15 | -0.119 | 0.211 | 0.496 |
| chelsa_bioclims_NZ.16 | 0.176 | 0.346 | -0.012 |
| chelsa_bioclims_NZ.17 | -0.234 | -0.243 | 0.252 |
| chelsa_bioclims_NZ.18 | -0.171 | 0.256 | 0.192 |
| chelsa_bioclims_NZ.19 | 0.116 | 0.151 | 0.289 |
Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
#get background env't for the species
nive_bg_env <- bg_env_crop(nive_loc,
species = "annulata",
environment = w,
buffer = 0.5)
#enfa for the sexual species
nive_sexual_enfa <- enfa_calc_fun(locs = nive_loc,
species = "annulata",
reproductive_mode = "sexual",
mask_raster = nive_bg_env)
#enfa for the asexual species
nive_asexual_enfa <- enfa_calc_fun(locs = nive_loc,
species = "annulata",
reproductive_mode = "asexual",
mask_raster = nive_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = nive_sexual_enfa$m,
asex_marg = nive_asexual_enfa$m,
full_species_name = "Niveaphasma annulata")###Spinotectarchus
summary.list.spin <- species_pca_fun(loc.clim, "spinotectarchus")
spin_plot <- plot_clim_pca(summary.list.spin$loc.clim, summary.list.spin$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(spin_plot, paste0(getwd(), "/plots/repro_mode_plots/spinotectarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
spin_loc <- loc %>%
filter(genus == "spinotectarchus")
#use sourced plot_locs_leaflet script to plot localities
spin_map <- plot_locs_leaflet(spin_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(spin_map, url = paste0(getwd(), "/plots/repro_mode_plots/spin_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/spin_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
Standard deviation 3.2696 2.2078 1.34284 0.86560 0.74380 0.4655 0.17853 0.16747 0.12847 0.11756
Proportion of Variance 0.5626 0.2566 0.09491 0.03944 0.02912 0.0114 0.00168 0.00148 0.00087 0.00073
Cumulative Proportion 0.5626 0.8192 0.91410 0.95354 0.98265 0.9941 0.99573 0.99721 0.99808 0.99881
PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.10182 0.07402 0.05291 0.04093 0.03100 0.02372 0.02118 0.01686 0.01123
Proportion of Variance 0.00055 0.00029 0.00015 0.00009 0.00005 0.00003 0.00002 0.00001 0.00001
Cumulative Proportion 0.99935 0.99964 0.99979 0.99987 0.99993 0.99995 0.99998 0.99999 1.00000
loadings.spin <- summary.list.spin$summary.pca$rotation
knitr::kable(round(loadings.spin[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | -0.279 | -0.110 | -0.223 |
| chelsa_bioclims_NZ.2 | -0.273 | -0.057 | -0.289 |
| chelsa_bioclims_NZ.3 | -0.277 | -0.148 | -0.178 |
| chelsa_bioclims_NZ.4 | 0.243 | -0.262 | -0.126 |
| chelsa_bioclims_NZ.5 | 0.199 | -0.325 | -0.112 |
| chelsa_bioclims_NZ.6 | 0.246 | -0.230 | -0.070 |
| chelsa_bioclims_NZ.7 | -0.133 | -0.273 | -0.032 |
| chelsa_bioclims_NZ.8 | 0.209 | -0.314 | -0.127 |
| chelsa_bioclims_NZ.9 | 0.254 | -0.217 | -0.106 |
| chelsa_bioclims_NZ.10 | 0.231 | -0.236 | -0.144 |
| chelsa_bioclims_NZ.11 | 0.215 | -0.303 | -0.135 |
| chelsa_bioclims_NZ.12 | 0.199 | 0.273 | -0.319 |
| chelsa_bioclims_NZ.13 | 0.174 | 0.242 | -0.434 |
| chelsa_bioclims_NZ.14 | 0.225 | 0.267 | -0.030 |
| chelsa_bioclims_NZ.15 | -0.165 | 0.138 | -0.579 |
| chelsa_bioclims_NZ.16 | -0.274 | -0.193 | -0.055 |
| chelsa_bioclims_NZ.17 | 0.217 | 0.281 | -0.234 |
| chelsa_bioclims_NZ.18 | -0.272 | -0.164 | -0.178 |
| chelsa_bioclims_NZ.19 | -0.206 | -0.046 | -0.140 |
###Tectarchus
summary.list.tect <- species_pca_fun(loc.clim, "tectarchus")
tect_plot <- plot_clim_pca(summary.list.tect$loc.clim, summary.list.tect$summary.pca, factor = "reproductive_mode")Ignoring unknown aesthetics: text
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. It's an annoying bug that hasn't been fixed yet.
#htmlwidgets::saveWidget(tect_plot, paste0(getwd(), "/plots/repro_mode_plots/tectarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
tect_loc <- loc %>%
filter(genus == "tectarchus")
#use sourced plot_locs_leaflet script to plot localities
tect_map <- plot_locs_leaflet(tect_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
#in case I want to save the map somewhere
#mapview::mapshot(tect_map, url = paste0(getwd(), "/plots/repro_mode_plots/tect_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/tect_map.pdf"))Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
Standard deviation 2.8210 2.5950 1.6213 0.97360 0.59328 0.41536 0.34340 0.23982 0.11480 0.08517
Proportion of Variance 0.4188 0.3544 0.1383 0.04989 0.01853 0.00908 0.00621 0.00303 0.00069 0.00038
Cumulative Proportion 0.4188 0.7733 0.9116 0.96150 0.98002 0.98910 0.99531 0.99833 0.99903 0.99941
PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19
Standard deviation 0.06494 0.04977 0.03772 0.03413 0.02872 0.02484 0.01571 0.01133 0.01044
Proportion of Variance 0.00022 0.00013 0.00007 0.00006 0.00004 0.00003 0.00001 0.00001 0.00001
Cumulative Proportion 0.99963 0.99976 0.99984 0.99990 0.99994 0.99997 0.99999 0.99999 1.00000
loadings.tect <- summary.list.tect$summary.pca$rotation
knitr::kable(round(loadings.tect[,1:3],3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| chelsa_bioclims_NZ.1 | 0.333 | 0.029 | -0.197 |
| chelsa_bioclims_NZ.2 | 0.310 | 0.019 | -0.288 |
| chelsa_bioclims_NZ.3 | 0.345 | 0.030 | -0.114 |
| chelsa_bioclims_NZ.4 | 0.048 | -0.380 | 0.021 |
| chelsa_bioclims_NZ.5 | 0.070 | -0.367 | 0.089 |
| chelsa_bioclims_NZ.6 | 0.030 | -0.379 | -0.011 |
| chelsa_bioclims_NZ.7 | 0.102 | 0.106 | 0.275 |
| chelsa_bioclims_NZ.8 | 0.071 | -0.369 | 0.078 |
| chelsa_bioclims_NZ.9 | 0.027 | -0.379 | -0.010 |
| chelsa_bioclims_NZ.10 | 0.018 | -0.379 | -0.015 |
| chelsa_bioclims_NZ.11 | 0.082 | -0.352 | 0.081 |
| chelsa_bioclims_NZ.12 | -0.279 | -0.060 | -0.352 |
| chelsa_bioclims_NZ.13 | -0.252 | -0.058 | -0.362 |
| chelsa_bioclims_NZ.14 | -0.300 | -0.045 | -0.285 |
| chelsa_bioclims_NZ.15 | 0.215 | -0.011 | -0.481 |
| chelsa_bioclims_NZ.16 | 0.351 | 0.042 | -0.021 |
| chelsa_bioclims_NZ.17 | -0.285 | -0.062 | -0.343 |
| chelsa_bioclims_NZ.18 | 0.273 | 0.008 | -0.209 |
| chelsa_bioclims_NZ.19 | 0.293 | 0.046 | -0.197 |
Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
This is for Tectarchus ovobessus.
#get background env't for the species
tect_bg_env <- bg_env_crop(tect_loc,
species = "ovobessus",
environment = w,
buffer = 0.5)
#enfa for the sexual species
tect_sexual_enfa <- enfa_calc_fun(locs = tect_loc,
species = "ovobessus",
reproductive_mode = "sexual",
mask_raster = tect_bg_env)
#enfa for the asexual species
tect_asexual_enfa <- enfa_calc_fun(locs = tect_loc,
species = "ovobessus",
reproductive_mode = "asexual",
mask_raster = tect_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = tect_sexual_enfa$m,
asex_marg = tect_asexual_enfa$m,
full_species_name = "Tectarchus ovobessus")This is an enfa for Tectarchus huttoni.
#get background env't for the species
tect_bg_env <- bg_env_crop(tect_loc,
species = "huttoni",
environment = w,
buffer = 0.5)
#enfa for the sexual species
tect_sexual_enfa <- enfa_calc_fun(locs = tect_loc,
species = "huttoni",
reproductive_mode = "sexual",
mask_raster = tect_bg_env)
#enfa for the asexual species
tect_asexual_enfa <- enfa_calc_fun(locs = tect_loc,
species = "huttoni",
reproductive_mode = "asexual",
mask_raster = tect_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = tect_sexual_enfa$m,
asex_marg = tect_asexual_enfa$m,
full_species_name = "Tectarchus huttoni")###Tepakiphasma Nothing. Only one locality.